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A Radiation Solver for the National Combustion Code 

Peter M. Sockol 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

A methodology is given that converts an existing finite volume radiative transfer method that requires 
input of local absorption coefficients to one that can treat a mixture of combustion gases and compute the 
coefficients on the fly from the local mixture properties. The Full-spectrum ^-distribution method is used 
to transform the radiative transfer equation (RTF) to an alternate wave number variable, g. The 
coefficients in the transformed equation are calculated at discrete temperatures and participating species 
mole fractions that span the values of the problem for each value of g. These results are stored in a table 
and interpolation is used to find the coefficients at every cell in the field. Finally, the transformed RTE is 
solved for each g and Gaussian quadrature is used to find the radiant heat flux throughout the field. The 
present implementation is in an existing cartesian/cylindrical grid radiative transfer code and the local 
mixture properties are given by a solution of the National Combustion Code (NCC) on the same grid. 
Based on this work the intention is to apply this method to an existing unstructured grid radiation code 
which can then be coupled directly to NCC. 


Introduction 

The National Combustion Code (NCC) is a state of the art Computational Fluid Dynamics (CFD) 
program specifically designed for combustion processes (Ref 1). The code employs an unstructured grid, 
massively parallel computing (Refs. 2 and 3), a dynamic wall function with the effect of adverse pressure 
gradient (Ref. 4), a low Reynolds number wall treatment (Ref 5), a cubic non-linear k-epsilon turbulence 
model (Ref 6), a lagrangian liquid phase spray model (Ref 7), and stiff laminar chemistry integration. 
Recently, viscous low-speed preconditioning (Ref 8) has been added to improve the low-speed 
convergence of the NCC in viscous regions and the ability to handle multiple sets of periodic boundary 
conditions has also been added. The combination of these features is usually not available in other CFD 
codes and gives the NCC an advantage when computing recirculating, turbulent, reacting spray flows. 
Previously, the NCC has undergone extensive validation studies for simple flows (Ref 9), complex flows 
(Ref 10), NOx emissions prediction (Ref 1 1) and traditional gas turbine combustor injectors (Ref. 12). 

For most high-temperature combustion problems, thermal radiation plays an important role in the heat 
transfer to the surrounding surfaces and the absorption and emission of radiation can also have a 
significant effect on the temperature distribution throughout the domain. Flence, if a designer wanted to 
use NCC to aid in the design of a jet engine combustor or a furnace combustion chamber, it would be 
very useful to be able to compute the radiative transfer as part of the process. In this work we have taken a 
significant step in that direction. 

At present two codes exist for solving the radiative transfer equation (RTE). The first is a stand-alone 
cartesian or cylindrical grid code (Ref 13) that solves in ID, 2D, or 3D. The second is an unstructured 
grid code (Ref 14) that could be coupled with NCC. Both of these codes require input of the local 
absorption coefficients. Hence the user of either code needs a means of coming up with at least 
approximate values of the absorption coefficients as a function of gas properties and wave number. 

A gas emits and absorbs radiation only at wave numbers in the vicinity of spectral lines. Since the 
main participating gases (CO 2 and H 2 O) each have more than one million absorption/emission lines, it 
would be impractical to solve the RTE at enough wave numbers to resolve this physics. Modest and 
coworkers have developed a methodology for handling this behavior by transforming the RTE to an 
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alternate wave number variable space where the coefficients vary smoothly (Ref 15). We currently have 
a Spectral Radiation Calculation Software package (Ref 16) (SRCS) which includes source code and data 
sets necessary to run the method as part of another application. 

In this work we use the existing cartesian/cylindrical radiation code modified to solve in the 
transformed wave number space. We also use the Modest full-spectrum k-distribution method and adapt 
the SRCS package to run with the modified cartesian/cylindrical radiation code. NCC is run on a cartesian 
grid to obtain local temperatures and species mole fractions. SRCS is used to obtain a table of absorption 
coefficients at discrete values of temperature and participating species mole fractions where these 
properties span the range of values from the NCC solution. We next interpolate within the table to obtain 
cell center values of k for every cell in the grid. Finally, we run the modified cartesian/cylindrical 
radiation code with these k values to obtain radiation intensity and heat transfer throughout the field. The 
expectation is that this same methodology will work quite well with the unstructured grid radiation code 
and thus provide a solver that can be coupled with NCC. 


Full-Spectrum A:-Distribution Method 

Over the years Modest and coworkers have added significant refinements to the basic method 
(Refs. 17 to 21). These increase the accuracy for inhomogeneous gas mixtures, but they also increase the 
complexity of the method. Since the turbulence- chemistry modeling used in NCC is still undergoing 
significant development, it was decided to use the basic full-spectrum k-distribution (FSK) method in this 
work (Ref. 15). In addition, since radiation scattering is not considered important for combustion gases, 
this effect is neglected (Ref 22). Also, while the presence of soot affects emission and absorption, 
scattering by soot is typically negligible. 

We begin with the spectral radiative transfer equation for an absorbing/emitting medium: 

( 1 ) 

where 7^ is the spectral intensity varying along a path s, 17 is the wavenumber, is the spectral 
absorption coefficient, 7^^ is the Planck function and 0 = {T, p, x) is an array of state variables. Here x is 
the mole fraction vector. We reorder the absorption coefficient in Equation (1) by multiplying by the delta 
function s{k — {rj, 0o) and integrating over the entire spectrum. Here /c^ {rf, is the absorption 
coefficient evaluated at some reference state 0o = (’^o.PoWo)- This gives 


^ = k* ( 0 , k) [/ {t. 00, k) 7, (T) - 7;,] , (2 

provided that at every wavenumber across the spectrum where /c^ ^ 0 q, p) has one and the same value 
k, we also have a unique value for /c^ ^0, 17 ^ = k* (^(p, k^. In Equation (2) 7^ and / are defined as 

4 = J IriS (k - Kr^ (00, 17 )) dp, (3 

/ (t, 00 , ^) = ]- J Ibv(T)S\k- Kr^ ( 00 , p)jdp. {4 
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In terms of the cumulative A:-distribution, 


rk r k 

g(T,^o,k)=l f(T,^o,f^)dk=l f (t , k*) dk* = g (t , k*Y 
*'0 *'0 


( 5 ) 


= k*(To,^,g'^ [a(T,To,g)Ii,(T) - Ig], ( 6 ) 

where 

Ig = h/f{To,±o,k), (7) 

a(T, To, g)=f (t, 0q, k)/ f (Tq, ^o,k), (8) 

and k* (j'q, <p, g^ is the k vs. g distribution as given by Equation (5), with the absorption coefficient 
evaluated at the local conditions 0 and the Planck function evaluated at the reference temperature Tq. 

The reordered spectral intensity 7^ can be obtained by solving the modified RTE, Equation (6), and the 
spectrally integrated intensity is then evaluated as 

00 pOO 

Igdg=l lkf(To,(l)o,k)dk=l Igdg. (9) 

I Jo Jo 



the RTE becomes: 



ds 


Since Ig is a smooth function of g, the integral in Equation (9) can be evaluated by Gaussian quadrature 
with values of /^computed and stored at the quadrature points. Similarly the integrated heat flux in any 
direction can also be found by Gaussian quadrature. 

Use of the SRCS package 

The SRCS package (Ref. 16) can be used to calculate the k versus g distribution as a function of the 
local state variables 0. However, to do this at every cell in a large 3D computation would be prohibitively 

expensive. Instead we construct a multidimensional table using discrete values of temperature and 
participating species mole fractions that span the values from the NCC solution. The SRCS software is 
used to populate the table with values of kg versus g^ where the g^ are the Gaussian quadrature points. 
We then interpolate within the table to obtain cell center values of kg for each gg at every spatial cell in 
the field. Note the same process is used to obtain Ug for each gg at every cell. 

We begin by calculating the reference state 0 q. Following Modest (Ref 15) we choose a volume 
averaged mole-fraction and a Planck-mean temperature: 

^ = y j xdV, (10) 
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( 11 ) 


Kp{To,^ = ^j Kp{T,x)li,(T)dV, 

and it is assumed that the pressure variation aeross the domain is negligible. For eombustor 
simulations, the volume V for these averages would be restrieted to regions with signifieant 
eoneentrations of eombustion produets (i.e., regions downstream of fuel injeetors). 

Next, for each (T, x ) grouping, we obtain the k (Tq.^Pi q) distribution using the SRCS software, 
interpolate within the distribution to obtain for each and store the kq in the table. Also, for each T 
in the table, we obtain a(T, Tq, g), interpolate to get for each g^ and store the in a second table. 
Since the absorption coefficients are expected to be more sensitive to T than to x, for the first table we 
construct cubic spline fits for kq versus T for each x and gq. We also construct cubic spline fits for 
versus T for each gq. 

Once this preliminary work has been completed, we proceed to compute kq and Uq for each cell in 
the field by interpolating within the two tables. Since for the present work we consider only two 
participating species (CO 2 and H 2 O), the spline fit is used to find kq for the local T at the four comers of a 
box containing the local x. Bilinear interpolation is then used to find the local value of kq. For three or 
four participating species this is readily extended to trilinear or quadralinear interpolation. We use the 
spline fit in the second table to find at the local T. 

Solution of the Radiative Transfer Equation 

The cartesian/cylindrical solver used in this study (Ref 13) uses a finite volume method to solve the 
transformed RTF, Equation (6). First note that Ig depends on both the spatial coordinates and the 
radiation direction (polar angle 9 and azimuthal angle 0). We integrate Equation (6) over a control 
volume AV and a control angle Afl and apply the divergence theorem to obtain 

[ [ Ig(s-n)dAdn= f f k*fTo,^,g)la(T,To,g)Ih(T)-Ig]dVdn (12) 

where s is the unit vector in the radiation direction and n is the outward unit normal along the control 
surface A A . If the intensity is assumed constant over the control volume AVi , then Equation (12) can be 
rewritten as 


Y^lljD\jAAtj= -B\l\+ S\. (13) 

J 

where the summation is over the faces of AVi , l\j is the intensity at face AAij of AV^ in direction , 


D-y = j (^^) 

/m' 

Bl = k-AViAn\ (15) 

Si = klaiOb)iAViAn\ (16) 
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and the subscript g has been dropped to simplify the notation. The face intensities llj are obtained from 
the step scheme. For D\j > 0 (an outflow of radiant energy), we use /• the center value of AVi . For D\j < 

0 (an inflow of radiant energy) we use /•, the center value of the adjacent control volume AVi, except 
where AAij lies on a domain boundary and a boundary condition is used to obtain ijj. Finally, the 
directional space is defined by taking constant spacing for 6 and 0, where 9 is measured relative to and 
(p is measured relative to in the — Cy plane. 

For a cartesian grid the solution of Equation (13) is straightforward. Start at one of the eight comers 
of the 3D domain and select that octant of the total directional space for which all the directions point into 
the interior of the grid. For each of these directions we sweep through the grid in such a way that all of the 
radiant inflow to a control volume comes either from adjacent control volumes that have been updated 
during the current sweep or from the adjacent boundary condition. Next move to another comer and 
repeat the process for a different octant of the directional space. After all eight comers have been used, 
the radiation intensity /• will have been updated for every control volume AV, and control angle AD.^ in 
the entire domain. 


Sample Cases 

Hot Gas Cylinder 

Before we present results from a combustion problem, we consider a simple case to show the kinds of 
results that can be obtained from the current cartesian/cylindrical radiation code. We choose a cylinder of 
hot gas with cold walls. The gas temperature Tg is set to 1000 K, the pressure ^ to 1 atm. and the mole 
fractions x^g of the participating gases to 0.2 and 0.02 for CO 2 and FI 2 O, respectively. The walls are set to 
0 K. Symmetry conditions are used to reduce the domain to radial distance r = 0 to rmax and axial 
distance z = 0 to 4 rmax with symmetry at r = 0 and z = 0. Even though the problem is axisymmetric, the 
azimuthal distance goes from = 0 to 2n. The number of spatial grid cells is taken as 20^ X 40r X lOOz 
with hyperbolic tangent stretching used to cluster the cells near the walls in r and z. The grid spacing in (p 
is kept uniform. The directional discretization uses 80 X 200 with uniform spacing in both directions. 
Finally we note that we solve the RTE for each of 10 Gaussian quadrature points Qg. The final 
distributions ofqy,qz,G and div(jq') are then found by Gaussian quadrature. 

For a first calculation we set rmax = 1.0 m. The resulting distributions are shown in Figure 1, where 
each of the variables is divided by an appropriate reference value: qref = oTg, Gref = 
qref/n, dvref = qref /rmax and a is the Stefan-Boltzman constant. In Figure 1(a) the radial heat 
flux rises gradually from zero at the centerline to a maximum at r = rmax, in Figure 1(b) the axial heat 
flux is close to zero until the end wall is approached at z = 4 rmax, in Fig. 1(c) the mean radiative 
intensity is largest at the centerline and decreases gradually as the cold walls are approached and finally in 
Figure 1(d) the radiative heat source is negligible until very near the cold walls. Thus this is an optically 
thick case. 

As a second calculation we set rmax = 0.01 m. The resulting distributions are shown in Figure 2, 
where the reference values are the same as in Figure 1 except that dvref is 100 times as large. The results 
shown in Figure 2(a) to (c) are similar to those in Figure 1. However, Figure 2(d) shows that the radiative 
heat source, while negligible near most of the centerline, gradually becomes substantial as the cold walls 
are approached. Thus this case could be considered optically thin. 

Sandia Flame D (Ref, 23) 

This is one of a series of piloted CHVair flames that were set up to provide data to test combustion 
models. In this case a mixture of 25% CH 4 and 75% dry air enters at a temperature of 294 K and mean 
velocity of 49.6 m/sec through the main jet located on the centerline with a diameter d = 12 mm. The jet 
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exit Reynolds number is 22400. The annular pilot, which surrounds the jet, has an outer diameter of 
18.2 mm. Burnt gas from the pilot enters at a mean temperature of 1880 K and a mean velocity of 
11.4 m/sec. Outside the pilot, dry air enters at a temperature of 29 1 K and a mean velocity of 0.9 m/sec. 
The entering composition of the pilot gas and the velocity and turbulence profdes, which were measured 
separately (Ref 24), are given in the above reference. 

This flame was simulated with NCC on a domain with radial distance r = 0 to 0.045 m and axial 
distance z = 0 to 0.375 m. A 2D-axisymmetric grid is used with 120 cells in r clustered at the outer edge 
of the jet, and 100 cells in z gradually stretching from the inlet to the exit. Axisymmetric flow was 
assumed. The calculations used a /c-epsilon turbulence model with only linear terms retained in the 
Reynolds stress. The chemistry model used finite rate chemistry with 12 species and 10 reactions and no 
turbulence-chemistry interaction. The resulting distributions are shown in Figure 3. Figure 3(a) shows hot 
temperatures reaching the centerline very close to the jet exit in contrast to the experimental data 
(Ref 23) where large temperatures along the centerline occur atz/c? >15. Figure 3(b) and (c) show that 
C02and FI 2 O mass fractions also reach large values on the centerline very close to the jet exit. The poor 
prediction for the flame location is likely a result of the lack of turbulence-chemistry interaction in these 
calculations. Since the purpose of this work is to demonstrate the capability of the radiation solver with 
input from NCC, this combustion simulation is still relevant. 

For the cartesian/cylindrical radiation solver, r and z are taken from the NCC grid and, although the 
problem is axisymmetric, the azimuthal distance cp again goes from ^ = 0 to 2n. The radiation solver 
spatial grid is taken as 20^ X 120r X lOOz and the directional grid is taken as 80 X 200 control angles. 
The local transformed absorption coefficents k* are calculated from the local values of temperature and 
mole fractions of the two participating species (CO 2 and FI 2 O) as obtained from the NCC solution. As in 
the above case we again solve the RTF for each of 1 0 Gaussian quadrature points and obtain the final 
distributions of , G and divfjq') by Gaussian quadrature. The resulting distributions are shown in 
Figure 4. Note the reference values are as defined for the Flot Gas Cylinder above with Tg set to 2000 K 
and rmax replaced by the jet diameter d = 0.0072 m. Since there are no solid walls, the temperatures on 
the domain boundaries are taken from the values at the adjacent grid cell and the emissivity on these 
boundaries is taken as 1.0. As a result Figure 4(a) and (b) shows that the radial and axial heat fluxes are 
relatively uninteresting. Figure 4(c) shows the radiant intensity increasing with z as the temperature 
increases. Figure 4(d) shows the radiative heat source largest in the interior of the flame along the 
centerline. 

Finally it should be noted that for this case the computer cost for the two codes is comparable. In 
particular, on an HP Z400 workstation, the time per iteration per cell was 1.8x10“^ sec for NCC while for 
the cartesian/cylindrical radiation code this time was 1.5x10“^ sec per spatial cell for a single solution of 
the RTF. Since in general it takes far fewer iterations to converge the RTF than it does for NCC, it would 
be quite practical to couple the two codes by running the radiation code every 10 to 100 iterations of NCC 
and then updating the radiation source term in the energy equation of NCC. This would allow the 
temperature distribution to adjust to the presence of the radiation coming from the flame. 


Concluding Remarks 

This paper presents a methodology for converting an existing finite volume radiative heat transfer 
method with input absorption coefficients to one that can treat a mixture of combustion gases and 
compute the coefficients on the fly from the local mixture properties. The Full- spectrum k-distribution 
method is used to transform the radiative transfer equation (RTF) to an alternate wave number variable, 
g. A Spectral Radiation Calculation Software package is used to calculate the transformed absorption 
coefficients at discrete values of temperatures and participating species mole fractions that span the values 
of the problem for each value of g. These results are stored in a table and interpolation is used to obtain 
cell centered values of the coefficients for every cell in the field. Finally, the transformed RTF is solved 
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for each g and Gaussian quadrature is used to find the radiant heat transfer, mean radiant intensity and 
radiant heat source throughout the field. 

In the present work the methodology is implemented in an existing cartesian/cylindrical grid radiation 
code and local mixture properties are obtained from a solution of the NCC on the same grid. The method 
is readily extendable to unstructured grids and provides the foundation for coupling radiation heat transfer 
with the flow solver in NCC. 
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Figure 1 . — Hot gas cylinder with cold walls: (a) radial heat flux, (b) axial heat flux, (c) 
mean radiant intensity and (d) heat flux divergence, optically thick case, rmax = 1.0m. 
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Figure 2. — Hot gas cylinder with cold walls: (a) radial heat flux, (b) axial heat flux, 
(c) mean radiant intensity and (d) heat flux divergence, optically thin case, 
rmax — 0.01m. 
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Figure 3. — Sandia Flame D — NCC results: (a) temperature distribution, (b) CO 2 mass 
fraction distribution, (c) FI 2 O mass fraction distribution, (d) CFI 4 mass fraction distribution. 
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Figure 4. — Sandia Flame D — radiation results: (a) radial heat flux, (b) axial heat flux, (c) mean radiant 
intensity and (d) heat flux divergence. 
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